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ABSTRACT 



We present a time-dependent multi-zone code for simulating the variability 
of Synchrotron- Self Compton (SSC) sources. The code adopts a multi-zone pipe 
geometry for the emission region, appropriate for simulating emission from a 
standing or propagating shock in a collimated jet. Variations in the injection of 
relativistic electrons in the inlet propagate along the length of the pipe cooling 
radiatively. Our code for the first time takes into account the non-local, time- 
retarded nature of synchrotron self-Compton (SSC) losses that are thought to 
be dominant in TeV blazars. The observed synchrotron and SSC emission is fol- 
lowed self-consistently taking into account light travel time delays. At any given 
time, the emitting portion of the pipe depends on the frequency and the nature 
of the variation followed. Our simulation employs only one additional physical 
parameter relative to one-zone models, that of the pipe length and is computa- 
tionally very efficient, using simplified expressions for the SSC processes. The 
code will be useful for observers modeling GLAST, TeV, and X-ray observations 
of SSC blazars. 

Subject headings: galaxies: active — quasars: general — radiation mechanisms: 
nonthermal — X-rays: galaxies 



1. Introduction 



In b lazars, radio loud active galaxies with their relativistic jets pointing close to our line 
of sig ht jBlandfordlEoTsI ). the observed radiation is dominated by relativistically beamed 
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emission from the sub-pc base of the jet. The blazar spectral energy distribution (SED) 
consists of two components. The first one, peaking at sub-mm to X-ray energies is al- 
most certainly due to synchrotron radiation, while the second one peaking at MeV to 
TeV energies is believed to be of inverse Compton (IC) nature, with both components pro- 
duced by the same population of relativistic electrons. The nature of the IC-scattered seed 
photo ns is still not clear, with both e xternal optical-UV photons from the broad line re- 
gion (ISikora. Begelman. fc Reeslll994l ) and IR photons from the putative molecular torus 
( iBlazejowski et al. II2OOOI ). as well as synchrotron photons (SSC, e.g. Maraschi, Ghisellini 
& Celotti 1992) contributing. It is believed that in the case of powerful blazars peaking at 
MeV to GeV energies, external seed photons from the broad line region dominate the IC 
scattering, while for weaker lineless blazars, peaking at ~ TeV energies, SSC is the dominant 
emission mechanism. Recent observational results (e.g. D'Arcangelo et al. 2007; Marscher 
et al. 2008), however, place the blazar emission site beyond the broad line region, lending 
support to the possibility that even in powerful blazars the GeV emission process may be 
pure SSC. For a review o f leptonic r nodels , as well as hadronic models for blazar emission 
(e.g. Aharonian 2000) see iBottcherl koO'^ . 



Due to the small angular size of the blazar emission region, it is not possible to spatially 
resolve the emitting region. Because of this, information about the structure of the emit- 
ting source can be obtained only through multiwavelength variability studies. Particularly 
telling is the variability of the emission produced by the highest energy electrons, because 
these electrons lose energy very quickly and exist only close to the sites where they have been 
produced. The goal of multiwavelength variability campaigns, involving in many cases ob- 
servations from radio up to TeV energies, is to study the characteristics of blazar variability, 
such as correlations and/or time delays between different energies, spectral characteristics of 
the observed variability, and the amplitude of variability as a function of energy. 

Most notable amongst blazars are the so called TeV blazars for which the synchrotron 
emission peaks at X-ray energies and the SSC emission peaks at TeV energies, as they 
present the active galaxies producing the highest confirmed electron energies. Variations 
of TeV blazars in these two bands can be extremely rapid (TeV doubling times as short 
as a few min; Aharonian et al. 2007), suggesting highly relativistic sub-pc scale flows 
(Doppler factors 6 ^ 50; e.g. Begelman, Fabian fc Rees 2008) that decelerat e substantially 
( Georganopoulos fc Kazanas 2003 : Ghisellini. Tayecchio fc Chiaberge 2005h to match the 
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much slower speeds required by VLBI observations (iPiner fc Edwards ll2004l : lPiner. Pant. Edwards 
20081 ). The TeV and X-ray variations are usually well correlated (e.g. Fossati et al. 2008; 
Maraschi et al. 1999; Sambruna et al. 2000), as expected, because they present variations 
by the same electron population. Usually, the lower energy emission within each of these 
bands peaks with small time delays relative to the higher energy emission (e.g. Fossati et al. 
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2000), while the X-ray and TeV spectra become harder with increasing flux (e.g. Takahashi 
et al. 1996). In certain cases, however, the X-ray and TeV variability do not seem to be 
correlated in a simple way (e.g. Aharonian et al. 2005). An intriguing variability pattern 
is that of the so-called 'orphan' flares, rare TeV flares that are not accompanied by X-ray 
flares (e.g. Krawczynski et al. 2004, Blazejowski et. al. 2005). While the correlated X-ray 
- TeV flares can be understood through an increase of the high energy emitting electrons, 
orphan TeV flares defy such a straightforward explanation. 

Models of blazar emission to date have, for the most part, been in some form of homo- 
geneous one zone models (e.g. Mastichiadis & Kirk 1997; Krawczynski, Coppi, & Aharonian 
2002). Such models, although appropriate for modehng the steady-state emission of a source, 
cannot simulate variability faster than the zone light crossing time. The basic limitation of 
one zone models stems from the fact that the high energy variability of both the synchrotron 
and SSC components is produced by high energy electrons with cooling times shorter than 
the light crossing time. Even if we assume that a disturbance in the radiating plasma (e.g. 
a higher density) instantaneously propagates across the zone, the received radiation would 
be smeared out for timescales shorter than the light crossing time, due to light travel time 
delays from different parts of the source ( §2.1j also Chiaberge & Ghiselhni 1999), and no 
variability faster than the light crossing time would be observed. One, therefore, cannot use 
one zone models to infer the source structure from high energy variability. 

Inhomogeneous variability models of increasing degree of sophistication have attempted 
to overcome the problems of one-zone models. The basic idea is to overcome the unphysical 
instantaneous injection throughout the source by adopting a specific geometry for the plasma 
flow that includes an inlet for injecting the radiating plasma. Variations in the injected 
plasma propagate and produce variations in the emissivity. Calculations of the received 
emission that take into account the light-travel times that radiation from different parts of 
the source takes to reach the observer, produce light curves that, at least, do not violate 
causality. How physically realistic these light curves are depends on the approximations used 
and on the characteristics of the source to be modeled. For example, synchrotron and IC 
losses from photons external to the source are local processes in the sense that, at a given 
point in the flow, the energy loss rate only depends on the local magnetic field and external 
photon field energy density, and not on the photon production throughout the source. This is 
not the case with SSC losses, because synchrotron photons produced throughout the source 
at earlier times - to take into account the light travel time from one point of the source 
to another - contribute to the photon energy density responsible for the SSC losses and to 
the emissivity at a given point and time in the source. To properly model sources like TeV 
blazars, in which SSC losses are important or even dominant, these considerations have to 
be taken into account. 
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There have been a f ew attempts during; t he last fifteen years to take these spatial con- 
siderations into account. iGomez et al. I (119941 ) considered a conical jet with a constant bulk 
Lorentz factor fiow in which the electron plasma and the magnetic field undergo adiabatic 
evolution onl y and calculated the radio va r iability induced by a shock wave propagating 
along the jet. iGeorganopoulos fc Marscher I (Il998a| ]bl) studied a parabohc jet that hydrody- 
namically accelerates and focuses to a conical geometry, and by following the synchrotron 
energy losses of the emitting electrons reproduced the radio to X-ray light curves of the 
X-ray bright blazar PKS 2155-304. This resulted in a frequency-dependent source size, in 
agreement with the fact that the variability timescale of synchrotron radiation increases with 
decreasing frequency. It also reproduced the usually observed soft lags (variations at soft 
X-rays being preceded by variations at hard X-rays) and the counterclockwise X-ray fiux 
- X-ray spectral index loops (e.g. Takahashi 1996, Maraschi 1999, Kataoka 2000, Ravasio 
2004), both manifestations of radiative cooling dominating the energetics of the high energy 
electrons. 



Kirk. Rieger. fc Mastichiadis I (119981 ) developed a semi-analytical model in which low 



energy electrons are injected in a zone where they undergo acceleration and eventually es- 
cape. The acceleration zone is assumed to move with a certain velocity, leaving behind 
the freshly accelerated electrons that cool through synchrotron radiation. Variations in the 
injection rate of low energy electrons in the acceleration zone result in variations of the 
emissivity, which are integrated over the volume of the source, taking into account time 
delays, to produce the observed multifrequency synchrotron light curves. This model in- 
cludes a treatment of particle acceleration and it is able to reproduce the uncommon hard 
lags (variations at hard X-rays p receded by variations at soft X-ra ys) and clockwise X-ray 
fiux - X-ray spectral index loops (jZhang II2002I : iRavasio et al. II2004I ). both manifestations of 
electrons still accelerating, just before reaching the maximum electron Lorentz factor, where 
the acceleration and radiative loss timescales are comparable. It does not include, however, 
SSC considerations. 



Chiaberge fc Ghisellini I (119991 ) studied the synchrotron and SSC emission, from a ho- 



mogeneous one zone model in which they assumed an instantaneous plasma injection, but 
taking into account the time delays with which the external observer would observe the 
variability (a similar approach was also taken by Kataoka et al. 2000). They also studied 
a case similar to that of Kirk et al. (1998), but without treating particle acceleration, by 
splitting the source into smaller one zone models that evolved autonomously, in the sense 
that (i) the SSC emission inside any one of their single zones uses as seed photons only the 
synchrotron photons produced in that zone and (ii) the SSC energy losses in every zone are 
caused only by the synchrotron photons produced in that zone. This simplified approach is a 
good approximation for following the energetics of the electrons if the source is synchrotron 
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dominated (because the SSC losses, although inappropriately calculated, are negligible), but 
does not produce realistic SSC light curves, because it does not calculate the emission due 
to upscattering synchrotron photons produced in other parts of the source in retarded times. 



A significant improvement was introduced by ISokolov. Marscher. fc McHardy I (120041 ) 
who incorporated in the calculation of the SSC emission from a given location in an inhomo- 
geneous source the synchrotron photons produced throughout the source in retarded times. 
This produces accurate SSC light curves, provided that the SSC losses that were sti ll treated 



as a local process are negligible. In a follow up paper, ISokolov fc Marscher I (120051 ). consid- 
ered also external Compton photons from the broad line region and the molecular torus. The 
challenge for inhomogeneous multi zone models for sources such as the TeV emitting blazars 
is the calculation of the non-local, time-retarded SSC losses induced by photons produced 
in other parts of the source. 

Here, we present such an inhomogeneous model that, for the first time, takes into 
account the non-local, time-delayed source emission on the SSC losses. We assume that a 
power law of relativistic electrons is injected at the inlet of a pipe, and that the electrons fiow 
downstream and cool radiatively. Variations in the injected electron distribution propagate 
downstream and manifest themselves as frequency dependent variability. This allows us to 
model high energy multiwavelength variability in a self-consistent manner. In §2] we describe 
the one zone model, which we use as a building block for the multizone model, and we 
show that, by construction, one zone models cannot simulate variability produced by high 
energy electrons with radiative cooling time shorter than the electron light crossing time 
from the single zone. In ^we describe our multizone, pipe-geometry model with emphasis 
on the coupling between subsequent zones and on the calculation of the local photon field 
due to non-local, time-delayed emission throughout the source. This is followed in ^ by a 
comparison of the code with analytical results and a series of case studies. We conclude in §21 
with a discussion of additional considerations that can be used as starting points for future 
work. 



2. The One Zone Model 

We consider a homogeneous spherical source of radius R permeated by a magnetic field 
B of energy density B'^/{87t). Energetic electrons are injected into the region at a rate 
ql'-fjt), where 7 is the electron Lorentz factor and t the injection time. These electrons lose 
energy through synchrotron and inverse Compton radiation and eventually escape after a 
characteristic time, tesc, of the order of the light crossing time. The implementation we 
describe is applicable to sources that are optically thin both to synchrotron emission in the 
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frequency range under consideration and to 7-ray absorption due to pair-production. 

The kinetic equation that describes the time-evolution of the electron energy distribution 
(EED) n(7,t) is 

Here, 7 includes both the synchrotron losses % and the inverse Compton losses jic in the 
Thomson regime (e7 < 3/4) 



4cT^ 2 . Aar 2 



min[emax, 3/(47)] 



7 = 7s + 7/c, = ^^^^^^c^ / U{e,t)de (2) 

where U{e, t) is the photon field energy density, e is the photon energy in units of the 
electron rest energy rrieC^, and o",- is the Thomson cross section. The photon field U{e,t) 
includes not only the synchrotron produced photons, but all the photons produced in the 
source through IC scattering, including, therefore, all the higher order SSC emission. 



We calculate the synchrotron emission following Melrose I (119801 ) : 



L,(e., t) = 1-85^-^ / z'/'e-Ml, t)dl, ^ = ( 3 j ^s/B.^', (3) 

where q is the electron charge, = B/Bcru, and Bcru = (mlc^) / (eh) is the critical magnetic 
field, where the electron cyclotron energy equals its rest mass and strong field considerations 
become important (e.g. Harding & Lai 2006). We note, that, although a (5-function approach 
for calculating the synchrotron emissivity would be faster, it would misinterpret the spectra in 
cases of hard power - law injection q oc 7"^, p > 2, whose cooling is known to produce a pile- 
up at its high energy cutoff (e.g. Kardashev 1962). Also, a 5-function synchrotron emissivity 
would not produce the oc u^^^ spectrum at frequencies below the critical frequency of the 
lowest energy electrons. These lowe r energy photons can be important seed photons for 



producing hard SSC TeV emission as iKatarzyhski et al. I (120061 ) point out 



To obtain the SSC emission through a simple integration as in the synchrotron case, we 
employ the ^-function approximation in which seed photons of energy eo are IC scattered by 
electrons of Lorentz factor 7 to energy ejc = (4/3)eo7^, as long as the scattering takes place 
in the Thomson regime (769 < 3/4). If we consider seed photons in the energy range de^ being 
IC scattered by electrons with Lorentz factors in the range d'y, then the emitted IC power is 
'yic'meC^n{'y,t)d'y and it is spread over a final photon energy range deic = 4(eoe/c/3)^''^(i7, 
resulting in an IC specific luminosity per seed photon energy interval 

dLjcieict) = m,c^n{-i,t)^ic-p- ^(7 - (3e7c/4eo)'/')e(3/4 - 760), (4) 

aeic 
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where in this context 7/c = (4cr^/3mc)7^f/(eo, t)(ieo and Q{x) is the Heaviside step function. 
This is written as 



dLjc{ejC; 



Integrating over the available seed photon distribution we obtain 



Lici^ic, t) 



n(7, t)U{eo, t)e~'^'S{j - (3e,c/4eo)'/')f^eo. (6) 



The range of final photon energies ejc is (4/3)e^eed,mm7mm < ^/c < Imax- For e/c within 
this range the limits of the above integration are: 



^seed,min ^IC — g ^seed,min^' 

Sejc , _ 4 

— ^ lOr e/c ^ 77 ^seed,minl' 



for eic < 7n 



2 

max 



(7) 
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L 4e 



ic 



for e/c > 7 



Following IChang &: Cooper I (1l970l ) and IChiaberge &: Ghisellini I (1l999l ). we discretise 
the kinetic equation ([T]), using a grid of logarithmically spaced Lorentz factors, 7^-, j = 
0,1,2, jmax, and linearly spaced time indices, ti. The difference equation that describes 
this system is 



At 



X ^ 7 

A7 tesc 



(9) 



Note that this is an implicit scheme, in the sense that the calculation of requires 
knowledge not of only the previous timestep EED, but also nj+i^j+i, the next higher 7 grid 
point at the current time. It is due to the implicit nature of the numerical procedure that 
this scheme is stable for large timesteps. The difference equation can be written as a system 
of tridiagonal equations 

n^-i+i = arij^i - + cqjA+i, (10) 



A7 , At . 

-r — , ... /. XT- ' = a— 7j+i,i+i, c = aAt. 

A7 + At A7 /tesc - At 7j- i+i A7 
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This can be easily computed if we use the initial condition rij^ = Vj (start with no 
relativistic electrons in the system) and the boundary condition Uj^^^^t = Vt. The first 
condition implies that the initial photon field is also zero for all photon energies. The 
second condition is satisfied if we set 7j,„„^ > 'Jmax, because the electrons can only lose 
energy, and there is no way to move to higher energies, populating the jmax bin of the 
7-grid. The simulation proceeds in the following manner: given n(7,tj) and U{eo,ti) we 
first calculate n(7,ti+i). We then calculate the synchrotron luminosity, Ls{e,ti^i), and the 
inverse Compton luminosity, Ljc{^,ti+i). The specific photon energy density f/(eo,ti+i) for 
the next timestep is obtained by adding these two luminosities and dividing by AnR'^meC^. 



2.1. Problems of one zone models in reproducing high energy variability 

By construction, in the one zone model variations in the injection propagate instanta- 
neously throughout the source, because no spatial coordinate enters the description of the 
system. If a power - law EED, q oc 7~^, 7^™. < 7 < Jmax, is injected in the source, radia- 
tive cooling and electron escape will result in a broken power law EED ^(7) in the source, 
steepening from an electron index p to p + 1, above 76, the electron Lorentz factor for which 
the escape time equals the radiative loss time 

where U stands for the total photon and magnetic field energy density in the source. It 
is these electrons, with 7 > 7?, that produce the high energy synchrotron and IC emission. 
Because the electron escape time is of the order of the light crossing time, t^sc = ktic = kR/c, 
k ~ 1— few, electrons with Lorentz factor 7 > 7ic = kjb, have a cooling time shorter than 
the light crossing time. One therefore anticipates that even for an injection event lasting 
much less that ti^, the high energy variable emission produced by electrons with 7 > 7;^ will 
be smeared out by light travel time effects and would appear to last for ~ tic, even though 
in each point in the source it lasts a shorter time ~ tididl- 

To demonstrate this, a flaring state was simulated by an increase by a factor of 5 in 
the injection g(7,t) that lasted tmj = tic/ 10. The system was allowed to reach a steady 
state before the disturbance in the injection was introduced. The emitted luminosity as a 
function of frequency was followed in time allowing us to produce light curves. By not taking 
into account time delays, and wrongly assuming that at any given time the observer sees 
the entire source as being at a single physical state, electrons with tcooi < tic produce high 
energy synchrotron and SSC variations that last less than tic (upper panel of Figured]). 



To take time delays into account, one has to consider that if at a certain time t the 
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observer receives photons from the nearest (front) part of the source, shces further away 
from the observer will be seen as they were in retarded times t — r/c, where r is the distance 
of the slice from the front part of the source. To treat this, at each time the luminosity 
of the system was recorded for a number of time steps covering a time equal to the light 
crossing time of the region. The luminosity observed at any time is thus the sum of the 
luminosity from each of these time steps, as each one represents the light emitted by a slice 
sequentially further back from the observer. The resulting light curves, plotted in the lower 
panel of Figure [H show that when the size of the region, and thus the time taken by light 
to travel across it, is accounted for, the observed variability of high energy electrons with 
trnni < Ur. Is spread out over t he length of the light crossing time. This is similar to what 



Chiaberge fc Ghisellini I (119991 ) observed when they performed a similar test. 



Another serious problem stemming from the lack of spatial considerations in the one- 
zone model comes from the fact that the model by construction assumes that the photons 
produced in the source at a given time are instantaneously available as seed photons for IC 
scattering throughout the source. This unphysical assumption has serious implications on 
the calculation of the SSC emissivity and on the calculation of the SSC losses, which in turn 
affect the evolution of the EED in the source, and through this the entire spectrum and light 
curves. The first effect has been addressed by the inhomogeneous model of Sokolov et al. 
(2004) and Sokolov & Marscher (2005). We present now the first multi-zone simulation that 
incorporates the issue of light travel time effects on the SSC losses. 



3. The Multi Zone Model 

3.1. The flow geometry 

The simplest and least computationally intensive deviation from a homogeneous model 
that can address the issues discussed above is one in which plasma is injected into a pipe 
of radius R and length L and cools radiatively as it flows downstream before it escapes 
after traversing the pipe length. Physically, this resembles the situation of a standing or 
propagating shock, as seen at the frame of the shock front. The plasma flow velocity u and 
the magnetic field B are constant along the pipe, and the EED is assumed to have no lateral 
gradients along the cross section of the pipe. Relativistic plasma is injected at the base of 
the flow. The injection variation timescale in this geometry can be arbitrarily smaller than 
R/c without violating causality, because, in principle, a disturbance in the plasma flow can 
reach the entire cross-section of the inlet at a single instance. However, the discretization 
procedure we describe below limits the range of meaningful injection variability to timescales 
greater than the plasma flow time through a zone of the flow. Our goal is to calculate the 
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EED in the frame of the pipe as a function of time and distance z from the inlet of the pipe, 
and through this calculate the emission received by an observer located at an angle Q to the 
axis of the pipe. 



3.2. The discretization of the pipe 

The pipe is broken down lengthwise and all cells are of length Z, comparable to the pipe 
radius R. The length of the pipe L = Nl, where N is the number of cells. Each cell is then 
simulated by a one zone model. The electron injection at the first cell is (3'i(7,t), similar to 
that defined for the one zone model. In each time step, the electrons that are calculated to 
leave each cell, are injected into the next cell in line in the next timestep. The injection of 
electrons, therefore, in cell i at time tj is qiipf, tj) = nj_i(7, tj-i)/tesc and the kinetic equation 
for the i-th cell is 

dnd^ + ^[7n.(7,t,)] + = (13) 

^esc ^esc 

Because in a time ~ tesc the electron content of a cell is transferred to the next cell, t^sc 
is connected to the bulk flow velocity u through tgsc = l/u. We also use tesc as the time 
step of our simulation. This ensures that the actual distance a disturbance in the electron 
distribution travels in a time step is equal to the bulk velocity times the time-step size, by 
transferring in a time step the electron content of cell i to cell i + The shortest variability 
timescale that can be simulated by this configuration is the single cell escape time. This 
limits the highest energy electrons that can be followed accurately to that of Lorentz factor 
7b = SrUeCu/AarUl, where U is the photon plus magnetic field energy density in the first cell 
(see equation [T2l) . and through this the highest energies of synchrotron and IC variations 
that can be reproduced. The advantage of the pipe configuration relative to a homogeneous 
model of the same size is that the highest energy electron variability we can follow is not 
connected to the length of the entire pipe, but to the length of a single zone, a quantity 
that is times shorter. This results in the pipe being able to track variations faster by A^, 
following electrons more energetic by A^, and synchrotron and SSC fequencies higher by A^^, 
relative to a homogeneous model of size L. An early version of this approach was presented 
by Graff et al. (2007). 



3.3. The photon energy density 



To solve the kinetic equation for each cell, an expression for the photon energy density 
resulting from all other cells by taking into account light travel time delays is required. In 
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general, for a region S characterized by a time-dependent emission coefficient j{r',t', e), tfie 
photon energy density U{r,t,e) is calculated by integrating j{T',t',e)/c in retarded times 
over the volume of the region S. Setting r = 0, for a point of interest in S, yields, without 
loss of generality, 

1 rr'in) 

U{r = 0,t,e) = - j{r',t' = t-r'/c,e)dr'dfl. (14) 

c Jo 

For our geometry we express this through the following approximation. Consider a cell i 
centered at Zi being illuminated by a cell m centered at Zm- The solid angle subtended at Zi 
by the cell m is AQ nR'^/ {zi — Zm^. The photon energy density AC/ {zi, Zm, t, e) at (zj, t) 
due to photons produced at Zj at retarded time t' — t — \zi — Zm\/c is: 

f.rr( . ^ j{Zm,t' = t-\Zi- Zm\/C,e) 7lR^ 

AU {zi, Zm, t, e) = — — -^l. (15) 

C yZi Zfyif 

Making use of the fact that the volume of each cell is Vc — irR^l, and that the luminosity 
L(e) emitted from a cell is L(e) = 47rj(e)V^, we obtain 

ATT/ 4. \ -^i^mit — t ~ \Zi — Zm\/C, e) (-\n\ 

= Anc{z. - z^y ■ ^''^ 

A summation over all cells in the pipe results to the total photon energy density U {zi, tj, e^) 
at cell i, time tj — jtesc, and energy 

TT( + x_ L{zi,tj,ek) ^ L{zm,t' = jtesc-\i-m\^/c,ek) , . 

where the first term is the photon energy density due to cell i itself. Note that a calculation 
of the photon energy density requires keeping record of a data-cube of the SED emitted by 
each cell at each time for a number of time steps equal to the light crossing time of the entire 
region. 



3.4. Pipe orientation 

When calculating the observed luminosities, we must take into account the different 
distances that light must travel from each of the different cells in the pipe to the observer. 
This difference is a function of the angle Q formed between the pipe and the observer. If 
photons emitted from the first cell are received by the observer at a given time, photons from 
cell i that are received simultaneously were emitted earlier by il cos 9 /c. Beaming can be 
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easily included in this model, by assuming that the entire pipe is moving with a relativistic 
velocity along its axis. Then for a choice of bulk Lorentz factor F and orientation angle dobs 
in the observer's frame, one can calculate the Doppler factor 6 and the angle 6 in the frame 
of the pipe, and transform the arrival times by dividing by 6, the observed frequencies by 
multiplying by 6, and the observed fluxes by multiplying by 6^. 

4. Results 

4.1. Comparison with analytical results for the synchrotron dominated case. 

A simple analytical test can be performed in the case of a source in which the energy 
losses are dominated by synchrotron radiation. In this case, adopting a steady power law 
electron injection at the inlet and assuming an electron residence time kL/c in the source, 
the source integrated EED will reach after time t » kL/c a steady-state. This steady- 
state, source-integrated electron distribution is a broken power law with an electron index 
steepening by one at ■jt, = i/imeC^) / {AarUskL). This will produce a synchrotron spectrum 
with a spectral break of 1/2 at an energy e^, = -B*7fe- For the synchrotron dominated 
configuration presented in Figure [2], the numerical result is in good agreement with the 
analytical both for the EED and the SED employing a J-function synchrotron emissivity. 
Note that the SSC luminosity is much lower than the synchrotron one. Note also that while 
the analytical synchrotron emission stops at = -B*7min5 the numerical continues to lower 
frequencies, due to the fu oc v^/"^ lower energy tail of the synchrotron emissivity. 

To compare the variability of our code with analytic expectations, we initiate injection of 
a power low EED at t = at the inlet and follow the evolution of the system toward a steady 
state. We expect that the EED will reach a steady state at or before a time ~ kL/c = 2L/c, 
equal to the time it takes for an electron to transverse the length of the pipe and escape, 
or equivalently the time it takes to fill the pipe with electrons. The part of the EED with 
7 < 7b will reach a steady state at ~ kL/c, because the electrons responsible for this emission 
transverse the entire pipe without cooling appreciably. At higher energies, 7 > 7f,, the 
electron radiative cooling lifetime t = [Sniec) / (AarUB'j) is shorter than t = kL/c. Electrons, 
therefore, of progressively higher energy will be confined closer to the inlet, resulting to an 
energy-dependent size pipe. This, together with the orientation of the observer, determines 
the time it takes for the emission at a given energy to reach the steady-state. 

For an observation angle 6 = 7r/2, there will be no position-dependent delays, given that 
the light path from all parts of the pipe to the observer are equal. Note that if the source 
is moving relativistically with bulk Lorentz factor F, 6* = 7r/2 transforms to ^ = 1/F in the 
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observers' frame. At < we expect a practically achromatic increase of the luminosity, 
reaching a steady state at ~ kL/c because the electrons responsible for this emission have 
7 < 7f,. At higher energies > the emission comes from electrons with energy 7 > 7;,, and 
the time to reach steady-state is t = (3mec)/(4cr^f/B7) = {^^'mecBl^'^) / {AarUB^l^"^) , where we 
have used = -8^,7^. 

To verify that the model variability agrees with our analytical predictions, we select 
four synchrotron energies, marked through the four vertical lines in Figure [2J The lowest 
energy (dotted vertical line) comes entirely from the z/^/^ tail of the synchrotron emissivity, 
and it is heavily dominated by the lowest energy electrons with 7 = 7mm << lb that have 
no time to cool before they escape. The second lower energy (broken vertical line in Figure 
12]) at ejnin < e < is predominately due to electrons with •ymin < 1 < Ih that also escape 
before they cool appreciably. In the upper panel of Figure [3] we plot the model light curves 
of these two low synchrotron energies, using the same line styles. As expected, the two light 
curves are almost indistinguishable, both reaching a steady state at t ~ 2L/c (marked by a 
solid vertical line in the upper panel of Figure [3]). 

The light curves of two more energies, this time with > e^, are marked by the dot- 
dash and triple dot-dash lines in Figure [2] and are plotted in the upper panel of Figure [3] 
using the same line styles. The vertical lines with the same line styles in Figure [3] indicate 
the time at which the corresponding light curves are expected to reach a steady state. As 
can be seen, at these times the light curves are at ~ 80% of their steady-state level. This is 
mainly because electrons continue to radiate at a given energy even when their Lorentz 
factor drops below (es/-Bv,)^/^, as the exponential decay of the emissivity indicates (equation 
[3l e.g. the synchrotron emissivity of an electron with Lorentz factor 7 = (es/S^)^/^ at time to 
drops oc exp[— (t/to)^] from its peak emissivity, requiring t = 2tQ to drop by 98%). To verifly 
this, we plot in the middle panel of Figure [3] the same light curves, using the ^-function 
approximation for the synchrotron emissivity. As can be seen, the two high energy light 
curves approach the steady-state level significantly closer to the analytically expected time. 

To evaluate if the light travel effects are properly taken into account, we rotate the pipe 
in such a way that the inlet is closer to the observer {6 = 0), as is the case for a propagating 
shock, that is observed 'jet on'. In this case, we anticipate the low energy (e < e^) variability 
that requires to fill the entire pipe with plasma, will reach a steady state after an additional 
time L/c, because this is the additional length that variations from the end of the pipe have 
to travel to reach the observer. The time it takes, therefore, for the lower energy emission 
steady-state to be reached is {k + l)L/c. For higher energy variability (e > e^) that reaches a 
steady-state before time t < kL/c, the additional light travel time required is {tc/k)/c = t/k, 
where tc/k is the maximum distance from the inlet of the pipe that the contributes to the 
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energy under consideration. The time it takes, therefore, for the steady-state to be reached 
is t{l + 1/k). We demonstrate these considerations in the lower panel of Figure O As 
can be seen the time it takes for the low energies (e < e?,) to reach steady state is now 
{k + l)L/c = 3L/c, and the time it takes for the high energy light curves to reach a steady 
state increases by a factor ~ 1/2 as can be seen by comparing the lower and upper panels 
of Figure [3] (e.g. the dot-dash light curve reaches 95% of the steady state after t ~ 1.1 L/c 
for 9 = 7r/2 and after t ^ 1.65 L/c for 6 = 0). 

An analytical result regarding the relation of the synchrotron to the SSC emission, that 
we can test our model against, addresses the relative amplitude of synchrotron and SSC 
emission. For states that are synchrotron dominated, an increase of the electron injection 
normalization results to a linear increase of the synchrotron emissivity because the syn- 
chrotron emissivity is proportional to the number of available electrons, and to a quadratic 
increase of the SSC luminosity, because the SSC luminosity is proportional to the product 
of number of electrons times the synchrotron photon energy density, which scales with the 
number of electrons (e.g. Ghisellini, Maraschi, & Dondi 1996). This result holds as long as 
the increased injection lasts long enough to occupy the entire volume of the source (in our 
case tyar > kL/c), bringing the source to a new steady state. 

To verify that our code reproduces this behavior we start from the configuration of 
Figure [21 which we observe from an angle 6 = tt/2 and, after we let the system reach its 
steady state, we increase the injected electron luminosity by a factor of 2 for t = 3L/c, longer 
than 2L/c, the time to reach steady state. The middle panel of Figure H] shows the evolution 
of the SED as the flare grows, with the times denoting time since the increased injection 
started, while the bottom panel shows the evolution of the flare as the flare dies out, starting 
from the time the additional injection is switched off. The characteristic quicker response 
of the high energy electrons is apparent. The two solid curves represent the low and high 
steady-states. At the upper panel we plot the light curves of the four frequencies marked 
with vertical lines at the two bottom panels. They are selected to roughly correspond to 
optical. X-ray, GeV and TeV energies, assuming that beaming will increase their observed 
values by a Doppler factor 5 ~ 20 — 40. Note that the synchrotron emission doubles when it 
reaches its steady-state, while the SSC quadruples, in agreement with the analytical result. 



4.2. Case studies 

We present here three variability case studies, for which we use as a starting point the 
configuration described in Figure [2l increasing the electron luminosity to L^j = 5 x 10^" 
erg s~^ to produce a steady-state SED (solid line in the lower panel of Figure E]) that for a 
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beaming 5 ~ 20 — 40 resembles those produced by flaring TeV blazars. We first study the 
case of an increased electron injection that lasts a short fraction of the light crossing time. 
After the system reaches its steady state, we increase the injected electron luminosity Linj by 
a factor of 2 for 0.1 L/c. At the bottom panel of Figure 0, we plot the SED evolution, while 
at the upper panel we plot the hght curves that correspond to the four energies denoted by 
vertical lines at the bottom panel. We note that the flare peaks with no time-delays at all 
energies, and then decays with the higher energies of each component decaying first. We also 
note that both the synchrotron and SSC flare have a higher amplitude at higher energies (and 
therefore the spectrum hardens as the flux increases). Also, because additional electrons are 
injected at all energies, the entire SED responds to the increased injection. The maximum 
fractional increase of the emission increases with frequency for both the synchrotron and the 
SSC components. This is because the higher the energy of the electrons required to produce 
a given synchrotron or SSC emission, the shorter their lifetime in the source; an additional 
injection, therefore, for a fixed time (O.lL/c in our case) will increase by a higher factor the 
number of higher energy electrons. 

To show how the variability event would be seen if we had the ability to resolve the 
pipe, we plot in the four panels of Figure the luminosity profile along the pipe for the four 
energies we study. In all cases we normalize the luminosity to the steady-state luminosity 
of the first zone. The lower curve in all cases depicts the steady-state luminosity profile. As 
expected, the high energy synchrotron and SSC emission are confined close to the inlet, while 
the low energy emission of both components extends throughout the source. The decline of 
the steady-state low energy SSC emission along the pipe is due to the fact that this emission 
is a convolution of a range of electron energies, and the higher energy electrons are gradually 
becoming unavailable away from the inlet. On top of the steady-states, we plot the snapshot 
luminosity profiles at the times depicted in the lower panel of Figure [51 It can be seen how the 
pulse is propagating away from the inlet, gradually disappearing due to cooling at the high 
synchrotron and SSC energies. At the low synchrotron and SSC energies it can clearly be 
seen how the pulse is spreading out due to the escape from zone to zone (the 71,(7, t)/tesc term 
of the kinetic equation). Note also that while the amplitude of the low energy synchrotron 
variation is substantial (starting with an increase by a factor of 2 at the inlet), because the 
pulse lasts for only O.lL/c, while low energy synchrotron emission is produced by electrons 
accumulating for 2L/c, the increase of the total low energy synchrotron emission is small 
(about 5%) as expected, and as can be seen in the upper panel of Figure O 

Another possible variation in the injected electron distribution is an increase in the 
maximum electron Lorentz factor of the EED, something that can result from a temporary 
increase in the electron acceleration rate. In this case the normalization of the EED remains 
constant and the increase in the injected luminosity depends on the electron index and the 
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new value of "jmax- Using the same configuration as above, we increase the value of "jmax 
by a factor of 5 for the same short time O.lL/c (for an electron index of p = 1.8 and for 
7mm = 10'^ this corresponds to an increase of Lj„j by ~ 50%) and we follow the evolution 
of the SED. As can be seen in Figure [3, the event is mostly manifested at the high energy 
tails of both the synchrotron and SSC components and dies quickly, because the high energy 
electrons injected have very short lifetimes. No significant variations are seen away from 
the high energy tails of the synchrotron and SSC components, a notable difference from the 
previous case in which the normalization of all electrons was increased. 

Finally, we study a case in which after the steady state is reached, an additional popu- 
lation of relatively low energy electrons is injected for a short time. This may be a plausible 
situation if one considers that a pre-acceleration mechanism is required to accelerate elec- 
trons up to Lorentz factors 7 ~ Tmp/me, where F is the typical bulk Lorentz factor of the 
flow, to provide the electrons that can be picked up by the Fermi acceleration mechanism 
for acceleration to much higher Lorentz factors (e.g. Sikora et al. 2002). Variations in this 
pre-acceleration mechanism that are not propagated to Fermi acceleration may account for 
the variations in the low energy tail of the electron distribution. We simulate this scenario 
by injecting an additional low energy FED with the same luminosity as the steady injection, 

but with -fmax = 10^. 

As can be seen in Figure [8|, as soon as the injection starts, the additionally injected 
low energy electrons produce low energy synchrotron emission (dotted line in Figure [H]) that 
remains at a plateau during the time it takes for the extra injection to transverse the length 
of the pipe, because these low energy electrons do not cool strongly. In the beginning, these 
low energy synchrotron photons are produced close to the TeV energy electrons of the steady 
flow (due to radiative losses these high energy electrons are found only close to the injection 
inlet) and are upscattered by them to TeV energies, producing additional TeV emission which 
is manifested as the rising part of the TeV flare (triple dot-dash line). As the event evolves, 
the additional population of low energy electrons propagates downstream and, due to the 
increasing distance between the additionally produced synchrotron photons and the pipe inlet 
where the TeV electrons are found, the TeV flare quickly subsides. The behavior of the X-ray 
emission (broken line) is very interesting: the extra low energy synchrotron photons produced 
increase the photon density experienced by the synchrotron X-ray emitting electrons and 
this causes additional cooling, slightly reducing the X-ray synchrotron emission. As the 
production of these additional low energy synchrotron photons is displaced downstream, 
their effect at the neighborhood of the inlet, where the synchrotron emitting electrons are 
found, decreases and the X-ray flux returns to its steady state. The behavior of the lower 
energy Gamma-ray emission does not show the sharp decline of the high energy Gamma 
rays. Instead, their light curve shows a gradual decline. This flare is mostly produced 
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by the additional electrons upscattering a broad energy range of seed photons that has a 
gradually decreasing level with distance from the inlet energy. The model behavior of a TeV 
f lare not accompanied by an X-ray flare is remini scent of the so-called orphan TeV flares 
( iKrawczynski et al I |2004J : iBlazejowski et al. II2005I ). These are rare flaring states of blazars 
characterized by an increase in the TeV luminosity that is not accompanied by a similar 
increase in X-ray energies. 



5. Discussion and future work 

We presented a multizone code that for the first time takes into account the non-local, 
time-retarded nature of SSC losses. This code is currently the only multizone model that 
incorporates the non-local, time-delayed SSC losses, and as such is uniquely suitable for 
modeling the results of multiwavelength campaigns at radio, optical. X-ray and 7-ray en- 
ergies, with the additional constraints in the critical and unexplored for TeV blazars GeV 
GLAST regime. 

As we argued, the results of one zone codes for the critical high energy regime of both 
the synchrotron and SSC components are problematic, and should not be used to infer the 
physical conditions in the source through variability modeling. We described our multizone 
code, tested it successfully against known analytical results, and presented a small number of 
variability case studies. The case studies we presented, although based on the same underly- 
ing steady-state configuration, exhibited very different variability patterns. This means that 
detailed modeling of broadband SEDs and simultaneous multiwavelength variability can be 
used to infer what is actually the cause of a given observed variability pattern, providing 
reliable constraints on the particle acceleration taking place. Orphan flares can be repro- 
duced assuming an increase of the injection of the low energy electrons, but not assuming 
the injection of a very high energy electron population, as we also showed analytically. 

The fact that this plausible variation cannot produce orphan flares significantly narrows 
the parameter space for events that can produce such events, possibly in agreement with 
their observed scarcity. 

The code we described can run with a typical workstation in a reasonable time of at 
most a few minutes at a resolution of ~ 10 bins per decade of observing frequency, ~ 10 
bins per decade of electron energy, and ~ 50 zones. To achieve this we employed a pipe 
geometry, and adopted an energy conserving 5-function approximation for the SSC emissivity, 
as well as a step function approach to take into account the change from the Thomson to 
Klein-Nishina IC scattering cross-sections. Adopting these approximations is problematic for 
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situations where IC scattering of narrow photon distributions (e.g. hne emission from the 
broad hne region or even a blackbody spectrum characterized by a typical photon energy eg) 
is important. In this case the adoption of the step function cross section description would 
create a strong artificial feature on the EED localized at the transition from the Thomson to 
the Klein-Nishina regimes at 7 oc 1/eo, which would then propagate to the emitted spectra 
through the ^-function IC emissivity. For SSC systems, however, where the seed photons are 
spread over many decades in energy, the resulting spectra are good approximations of those 
produced using the full expressions for the synchrotron and SSC emissivities as well as the 
full Klein-Nishina cross section. 

Including the above considerations, as well as the processes of synchrotron opacity and 
pair production through 7-ray absorption within the source, would increase the execution 
time up to levels marginally comfortable for the typical workstation. Most probably, such 
an extension of the code would require parallelization. A more desirable upgrade of the 
code would drop the assumption of no lateral gradients in the plasma characteristics by 
switching to a two-dimensional geometry, in which the electron distribution and the SSC 
photon energy density are allowed to change laterally to the flow direction. Such consid- 
erations may be relevant to the r ecently observed ~ 0. 75 days delay between the IR and 



the X-ray variability in 3C 273 (IMcHardy et al. 1120071 ). These authors argued that the 



delay may be attributed to the time it takes for the SSC photon energy density to built 
up as the SSC photons are transversing the cross section of the flow. This upgrade will 
scale the computation time roughly by A^^, increasing it from ~ few minutes to ~ sev- 
eral hours. We note here tha t our formalism can be extended to treat velocity profiles in 
term of the decelerating flow jGeorganopoulos fc Kazanas 2003 ) or the spine sheath model 



(jGhisellini. Tavecchio fc Chiaberge II2005I ) that have developed to address the lack of super- 



luminal motions in TeV blazars (e.g. Finer & Edwards 2004). 

Another upgrade that can be incorporated in the existing code, this time with a minimal 
co mputational overhead, is that o f a zo ne for particle acceleration, following the formalism 
of iKirk. Rieger. fc Mastichiadis I Jl998h . In this case, in the first zone of the model, low 



energy electrons will be injected and allowed to accelerate while suffering radiative losses 
due to synchrotron and non-local SSC. These particles will subsequently escape into the 
pipe and flow downstream. This configuration will require a different numerical scheme 
for the acceleration zone, since there most particles are advected upward in energy space, 
but there is a possibility, in a time-dependent scenario, of the highest energy particles being 
advected downward, while the rest of the electrons are still advected upwards. The benefit of 
including particle acceleration in the code is that it will allow us to study cases of hard lags/ 
counterclockwise loops in the X-ray hardness - X-ray flux diagrams thought to result when 
acceleration and loss timescales are comparable. Such a code could be used to model the 



observed curved X-ra y spectra of high pea k frequency blazars in the framework of episodic 
particle acceleration (jPerlman et al. II2005I ). 
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Fig. 1. — A short variation by a factor of 5 for = tic/10, 1^. = R/c, is introduced in the 
one zone model. The dotted Une tracks the synchrotron emission of the low energy electrons 
with cooling time tcooi > tesc, the broken line of electrons with tic < icooi < tesc, the dash- 
dot line of electrons with tinj < tcooi < Uc, ^-nd the solid line of electrons with tcooi < Unj- 
In the upper frame, no light crossing delays are taken into account and this results to the 
unphysical result of variability times shorter than tic, for radiation produced by electrons 
with tcooi < tic- The lower frame depicts the same event with light crossing delays taken into 
account. In this case, the light crossing time is the smallest observable variability scale. The 
time-integrated emission in the disturbance is the same in both cases, as can be easily seen 
for the highest frequency that reaches a plateau in both panels: (5 — 1) x 0.1 = (1.4 — 1) x 1. 
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Fig. 2. — Comparison of the steady-steady numerical result with the analytical solution of 
a synchrotron dominated configuration. The following parameters have been used: A pipe 
of length L = 10^^ cm and radius R = 5 x 10^^ cm (aspect ratio 10 : 1). The magnetic field 
in the pipe is i? = 0.3 G, and a power law EED is continuously injected with •jrnin = 10^, 
Imax = 2 X 10®, p = 1.8. The pipe is split into 40 cylindrical slices of equal height / = L/40 = 
2.5 X lO^^cm, and the escape time is set to = 2 times the light crossing time. With these 
parameters fixed, the ratio of the SSC to synchrotron luminosity increases with increasing 
electron injected power, and for Linj = 5 x 10^^ erg s~^ the system is synchrotron dominated. 
At the lower panel we plot the analytical (broken line) and numerical (solid line) steady-state 
EED. At the upper panel we plot the analytical synchrotron (broken line) and numerical 
(solid line) synchrotron and SSC SED. The four vertical lines mark four frequencies for which 
we study their variability in Figure [3] using the same line styles. 
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Fig. 3. — The light curves of the four frequencies shown by the vertical lines in Figure [21 using 
the same line styles, for two different pipe orientations. The solid vertical lines correspond 
to the time kL/c required for the low energy (e < e^) light curves to reach steady state, 
while the dot-dash and triple dot-dash vertical lines correspond to the analytical estimate 
for the time it takes for the high energy (e > ei,) hght curves to reach steady state. The pipe 
orientation is = 7r/2 for the upper and middle panels, and = for the bottom panel. 
The middle panel uses a 5-function approximation for the synchrotron emissivity. 
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Fig. 4. — A quadratic variation: a doubling of the injected electron power for t = 3L/c, 
a time longer than 2L/c, the time it takes electrons to transverse the pipe. The initial 
configuration is the same as that of Figure O The middle panel shows snapshots of the SED 
for a range of times elapsed from the beginning of the additional injection and the bottom 
panel shows snapshots of the SED evolution after the additonal injection has been switched 
off. The upper panel shows the light curves for the four energies depicted by vertical lines 
at the lower two panels. 
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Fig. 5. — Here we use the same steady-state configuration as in Figure [2l except for fiiglier 
injected EED luminosity, Linj = 5 x 10'^° erg s~^. After tfie system readies tlie steady-state, 
we increase the injected luminosity by a factor of 2 for 0.1 L/c. The bottom panel shows the 
evolution of the SED, while the upper panel the light curves of the flare for the four energies 
marked by the vertical lines at the bottom panel. 
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Fig. 6. — The luminosity profile along the pipe for the four energies whose light curves are 
plotted in Figure 0, for the same pipe configuration. In each case the lower curve represents 
the luminosity profile, normalized to the steady-state luminosity of the first zone. The 
profiles at times t = 0.1, 0.2, 0.3, 0.5, 1.0 light crossing times are also plotted and can be 
seen moving away from the inlet as the additional injection propagates. 
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Fig. 7. — Starting from the same steady state as in Figure [21 •jmax is increased by a factor 
of 5 for O.lL/c. 
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Fig. 8. — Starting from the same steady state as in Figure[2l we inject for O.IL/ c an additional 
low energy EED with the same luminosity as the steady injection, but with 'jmax = 10"^. 



